DLPFC brain development

Setup

library('SummarizedExperiment')

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, cbind, colMeans,
##     colnames, colSums, do.call, duplicated, eval, evalq, Filter,
##     Find, get, grep, grepl, intersect, is.unsorted, lapply,
##     lengths, Map, mapply, match, mget, order, paste, pmax,
##     pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
##     rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which, which.max, which.min

## Loading required package: S4Vectors

## Warning: package 'S4Vectors' was built under R version 3.4.2

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## Loading required package: GenomeInfoDb

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: DelayedArray

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following object is masked from 'package:base':
## 
##     apply

library('recount')

## Warning: package 'recount' was built under R version 3.4.2

## Setting options('download.file.method.GEOquery'='auto')

## Setting options('GEOquery.inmemory.gpl'=FALSE)

library('edge')
library('devtools')
library('ggplot2')
library('reshape2')

load('../genotype_data/mds_extracted_from_BrainSeq_Phase2_RiboZero_Genotypes_n546.Rdata')

rse <- lapply(c('../count_data/dlpfc_ribozero_brainseq_phase2_hg38_rseGene_merged_n449.rda', '../count_data/hippo_brainseq_phase2_hg38_rseGene_merged_n442.rda'), function(f) {
    load(f)
    rowRanges(rse_gene)$meanExprs <- NA
    
    ## Add mds info
    m <- match(colData(rse_gene)$BrNum, rownames(mds))
    colData(rse_gene) <- cbind(colData(rse_gene), mds[m, ])
    return(rse_gene)
})

rse <- lapply(rse, function(r) { r[, colData(r)$Dx == 'Control']})
rpkm <- lapply(rse, recount::getRPKM, 'Length')

rpkm_all <- do.call(cbind, rpkm)

meanRPKM <- rowMeans(rpkm_all)
keepIdx <- meanRPKM > 0.5
table(keepIdx)

## keepIdx
## FALSE  TRUE 
## 37022 21015

colData(rse[[1]])$Region <- 'DLPFC'
colData(rse[[2]])$Region <- 'HIPPO'

## Define models
fm_mod <-  ~Age + fetal + birth + infant + child + teen + adult + Sex + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
fm_mod0 <- ~ Age + Sex + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
fm_mod_all <- ~Age *Region + fetal * Region + birth *Region + infant *Region + child * Region + teen * Region + adult * Region + Sex + Region + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
fm_mod0_all <- ~ Age + fetal + birth + infant + child + teen + adult + Sex + Region + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5


get_mods <- function(pd, int = FALSE, only_pd = FALSE) {
    
    if(int) {
        pd <- rbind(
            pd[[1]][, c('Age', 'Race', 'Sex', 'Region', 'snpPC1', 'snpPC2', 'snpPC3', 'snpPC4', 'snpPC5')],
            pd[[2]][, c('Age', 'Race', 'Sex', 'Region', 'snpPC1', 'snpPC2', 'snpPC3', 'snpPC4', 'snpPC5')])
        pd$Region <- relevel(factor(pd$Region), ref = 'DLPFC')
    }
    
    fetal = ifelse(pd$Age < 0, 1,0)
    birth = pd$Age
    birth[birth < 0] = 0 # linear spline
    infant = pd$Age - 1
    infant[infant < 0] = 0 # linear spline
    child = pd$Age - 10
    child[child < 0] = 0 # linear spline
    teen = pd$Age - 20
    teen[teen < 0] = 0 # linear spline
    adult = pd$Age - 50
    adult[adult < 0] = 0 # linear spline
    
    pd$Race <- relevel(factor(pd$Race), ref = 'CAUC')
    pd$Sex <- relevel(factor(pd$Sex), ref = 'F')
    
    pd$fetal <- fetal
    pd$birth <- birth
    pd$infant <- infant
    pd$child <- child
    pd$teen <- teen
    pd$adult <- adult
    
    if(only_pd) {
        return(as.data.frame(pd))
    }
    
    ### adjust for race, sex
    if(int) {
        mod = model.matrix(fm_mod_all, data=pd)
        mod0 = model.matrix(fm_mod0_all, data=pd)
    } else {
        mod = model.matrix(fm_mod, data=pd)
        mod0 = model.matrix(fm_mod0, data=pd)
    }
    
         
    return(list(mod = mod, mod0 = mod0))
}

DLPFC

### DLPFC

if(!file.exists('edge_dlpfc.Rdata')) {
    de_obj_dlpfc <- build_models(data = rpkm[[1]][keepIdx, ], full.model = fm_mod, null.model = fm_mod0, cov = get_mods(colData(rse[[1]]), only_pd = TRUE))
    ef_obj_dlpfc <- fit_models(de_obj_dlpfc, stat.type = "lrt")
    de_lrt_dlpfc <- lrt(de_obj_dlpfc, nullDistn = "normal")
    system.time( de_odp_dlpfc <- odp(de_obj_dlpfc, bs.its = 100, verbose = FALSE, n.mods = 50) )
    save(de_obj_dlpfc, ef_obj_dlpfc, de_lrt_dlpfc, de_odp_dlpfc, file = 'edge_dlpfc.Rdata')
} else {
    load('edge_dlpfc.Rdata')
}

summary(de_lrt_dlpfc)

## 
## ExpressionSet Summary 
##  
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 21015 features, 296 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: R10424 R12258 ... R5885 (296 total)
##   varLabels: SAMPLE_ID FQCbasicStats ... adult (83 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:  
## 
## de Analysis Summary 
##  
## Total number of arrays: 296 
## Total number of probes: 21015 
##  
## Biological variables: 
##  Null Model:~Age + Sex + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
##  Full Model:~Age + fetal + birth + infant + child + teen + adult + Sex + 
##     snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
## ....... 
##  
## 
## Statistical significance summary:
## pi0: 0.09010914  
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05  <0.1    <1
## p-value    15479  16617 17854  18371 18809 19248 20976
## q-value    16555  17811 19258  19927 20421 21015 21015
## local fdr  14833  15918 17159  17677 18095 18523 20134

hist(qvalueObj(de_lrt_dlpfc))

summary(qvalueObj(de_lrt_dlpfc)$qvalues)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.000e+00 0.000e+00 0.000e+00 3.711e-03 2.258e-05 9.011e-02

summary(de_odp_dlpfc)

## 
## ExpressionSet Summary 
##  
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 21015 features, 296 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: R10424 R12258 ... R5885 (296 total)
##   varLabels: SAMPLE_ID FQCbasicStats ... adult (83 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:  
## 
## de Analysis Summary 
##  
## Total number of arrays: 296 
## Total number of probes: 21015 
##  
## Biological variables: 
##  Null Model:~Age + Sex + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
##  Full Model:~Age + fetal + birth + infant + child + teen + adult + Sex + 
##     snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
## ....... 
##  
## 
## Statistical significance summary:
## pi0: 0.7455312   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05  <0.1    <1
## p-value     9798  11757 13903  14427 14794 15297 21015
## q-value     9738  11095 13835  14375 14758 15272 21015
## local fdr   9619   9694 10740  12610 13211 13537 14691

hist(qvalueObj(de_odp_dlpfc))

summary(qvalueObj(de_odp_dlpfc)$qvalues)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000008 0.0000008 0.0006612 0.1302452 0.1576178 0.7106748

HIPPO

### HIPPO

if(!file.exists('edge_hippo.Rdata')) {
    de_obj_hippo <- build_models(data = rpkm[[2]][keepIdx, ], full.model = fm_mod, null.model = fm_mod0, cov = get_mods(colData(rse[[2]]), only_pd = TRUE))
    ef_obj_hippo <- fit_models(de_obj_hippo, stat.type = "lrt")
    de_lrt_hippo <- lrt(de_obj_hippo, nullDistn = "normal")
    system.time( de_odp_hippo <- odp(de_obj_hippo, bs.its = 100, verbose = FALSE, n.mods = 50) )
    save(de_obj_hippo, ef_obj_hippo, de_lrt_hippo, de_odp_hippo, file = 'edge_hippo.Rdata')
} else {
    load('edge_hippo.Rdata')
}

summary(de_lrt_hippo)

## 
## ExpressionSet Summary 
##  
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 21015 features, 309 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: R11135 R11137 ... R5768 (309 total)
##   varLabels: SAMPLE_ID RNum ... adult (74 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:  
## 
## de Analysis Summary 
##  
## Total number of arrays: 309 
## Total number of probes: 21015 
##  
## Biological variables: 
##  Null Model:~Age + Sex + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
##  Full Model:~Age + fetal + birth + infant + child + teen + adult + Sex + 
##     snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
## ....... 
##  
## 
## Statistical significance summary:
## pi0: 0.028477    
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05  <0.1    <1
## p-value    16429  17622 18765  19214 19526 19840 21015
## q-value    18162  19336 20498  20943 21015 21015 21015
## local fdr  16307  17532 18677  19118 19427 19754 20873

hist(qvalueObj(de_lrt_hippo))

summary(qvalueObj(de_lrt_hippo)$qvalues)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.000e+00 0.000e+00 0.000e+00 6.873e-04 8.690e-07 2.847e-02

summary(de_odp_hippo)

## 
## ExpressionSet Summary 
##  
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 21015 features, 309 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: R11135 R11137 ... R5768 (309 total)
##   varLabels: SAMPLE_ID RNum ... adult (74 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:  
## 
## de Analysis Summary 
##  
## Total number of arrays: 309 
## Total number of probes: 21015 
##  
## Biological variables: 
##  Null Model:~Age + Sex + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
##  Full Model:~Age + fetal + birth + infant + child + teen + adult + Sex + 
##     snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
## ....... 
##  
## 
## Statistical significance summary:
## pi0: 0.3230261   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05  <0.1    <1
## p-value    15335  16496 17652  18127 18500 18831 21015
## q-value    15748  16991 18159  18659 19003 19334 21015
## local fdr  14290  15182 16358  16838 17175 17499 18693

hist(qvalueObj(de_odp_hippo))

summary(qvalueObj(de_odp_hippo)$qvalues)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 2.600e-07 2.600e-07 2.600e-07 2.438e-02 1.023e-04 3.055e-01

Quick region comparison

### Compare quickly

stopifnot(identical(rownames(de_obj_dlpfc), rownames(de_obj_hippo)))

comp <- data.frame(
    DLPFC = qvalueObj(de_lrt_dlpfc)$qvalues,
    HIPPO = qvalueObj(de_lrt_hippo)$qvalues,
    DLPFC_odp = qvalueObj(de_odp_dlpfc)$qvalues,
    HIPPO_odp = qvalueObj(de_odp_hippo)$qvalues,
    gene_id = rownames(de_obj_dlpfc),
    stringsAsFactors = FALSE
)
dim(comp)

## [1] 21015     5

## Using the LRT model results
table(comp$DLPFC < 1e-04, comp$HIPPO < 1e-04)

##        
##         FALSE  TRUE
##   FALSE  1228  3232
##   TRUE   1625 14930

## Using the ODP results
table(comp$DLPFC_odp < 1e-04, comp$HIPPO_odp < 1e-04)

##        
##         FALSE TRUE
##   FALSE  3966 7311
##   TRUE   1301 8437

save(comp, file = 'edge_comp.Rdata')

Both regions

### Both
mods <-  get_mods( list(colData(rse[[1]]), colData(rse[[2]])), int = TRUE)
sapply(mods, colnames)

## $mod
##  [1] "(Intercept)"        "Age"                "RegionHIPPO"       
##  [4] "fetal"              "birth"              "infant"            
##  [7] "child"              "teen"               "adult"             
## [10] "SexM"               "snpPC1"             "snpPC2"            
## [13] "snpPC3"             "snpPC4"             "snpPC5"            
## [16] "Age:RegionHIPPO"    "RegionHIPPO:fetal"  "RegionHIPPO:birth" 
## [19] "RegionHIPPO:infant" "RegionHIPPO:child"  "RegionHIPPO:teen"  
## [22] "RegionHIPPO:adult" 
## 
## $mod0
##  [1] "(Intercept)" "Age"         "fetal"       "birth"       "infant"     
##  [6] "child"       "teen"        "adult"       "SexM"        "RegionHIPPO"
## [11] "snpPC1"      "snpPC2"      "snpPC3"      "snpPC4"      "snpPC5"

pd_all <- get_mods( list(colData(rse[[1]]), colData(rse[[2]])), int = TRUE, only_pd = TRUE)

if(!file.exists('edge_both.Rdata')) {

    de_obj <- build_models(data = rpkm_all[keepIdx, ], full.model = fm_mod_all, null.model = fm_mod0_all, cov = pd_all)
    system.time( ef_obj <- fit_models(de_obj, stat.type = "lrt") )
    system.time( de_lrt <- lrt(de_obj, nullDistn = "normal") )
    
    ## odp.stat and null.stat are full of NaNs, which eventually
    ## leads to qvalue() to fail. 
    # de.fit <- fit_models(de_obj, stat.type = 'odp')
    # odp.parms <- kl_clust(de_obj, de.fit, n.mods = 50)
    # odp.stat <- edge:::odpStat(n.res = de.fit@res.null, clustParms = odp.parms)
    # null.stat <- edge:::bootstrap(object = de_obj, obs.fit = de.fit, clustParms = odp.parms, bs.its = 2, verbose = FALSE)
    # pval <- qvalue::empPvals(stat = odp.stat, stat0 = null.stat)
    # qval <- qvalue::qvalue(pval)
    #
    
    #system.time( de_odp <- odp(de_obj, bs.its = 100, verbose = FALSE, n.mods = 50) )
    
    #save(de_obj, ef_obj, de_lrt, de_odp, file = 'edge_both.Rdata')
    save(de_obj, ef_obj, de_lrt, file = 'edge_both.Rdata')
} else {
    load('edge_both.Rdata')
}

summary(de_lrt)

## 
## ExpressionSet Summary 
##  
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 21015 features, 605 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: R10424 R12258 ... R5768 (605 total)
##   varLabels: Age Race ... adult (15 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:  
## 
## de Analysis Summary 
##  
## Total number of arrays: 605 
## Total number of probes: 21015 
##  
## Biological variables: 
##  Null Model:~Age + fetal + birth + infant + child + teen + adult + Sex + 
##     Region + snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
##  Full Model:~Age * Region + fetal * Region + birth * Region + infant * Region + 
##     child * Region + teen * Region + adult * Region + Sex + Region + 
##     snpPC1 + snpPC2 + snpPC3 + snpPC4 + snpPC5
## 
## ....... 
##  
## 
## Statistical significance summary:
## pi0: 0.3320878   
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05  <0.1    <1
## p-value    11767  13387 15260  16081 16791 17552 20772
## q-value    12135  13951 16019  17004 17815 18671 21015
## local fdr  10312  11710 13424  14239 14815 15486 18098

hist(qvalueObj(de_lrt))

summary(qvalueObj(de_lrt)$qvalues)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000000 0.0000000 0.0000076 0.0307749 0.0077139 0.3320878

#summary(de_odp)
#hist(qvalueObj(de_odp))
#summary(qvalueObj(de_odp)$qvalues)

Explore a gene

The following code explores a single gene based on the de_lrt_dlpfc, de_odp_dlpfc, de_lrt_hippo, de_odp_hippo and de_lrt results.

gene_explore <- function(ef, de, pd) {
    fitVals <- fitFull(ef)
    fitVals0 <- fitNull(ef)

    i <- which.min(qvalueObj(de)$qvalue)
    print(qvalueObj(de)$qvalue[i])

    df <- data.frame(age= pd$Age, sex= pd$Sex, race = pd$Race, region = pd$Region, null.model = fitVals0[i, ], full.model = fitVals[i, ], raw = exprs(de)[i, ])
    
    
    df <- melt(data = df, id.vars=c("age", "sex", 'race', 'region'))
    df$value <- log(df$value + 1)
    
    p1 <- ggplot(df, aes(log10(age + 1), value, color=sex, shape = race)) + geom_point() + facet_wrap(region~variable) + ylab('log(RPKM + 1)') + xlab('log10(age + 1)') + geom_vline(xintercept = log10(c(0, 1, 10, 20, 50) + 1))
    p2 <- ggplot(df, aes(age, value, color=sex, shape = race)) + geom_point() + facet_wrap(region~variable) + ylab('log(RPKM + 1)') + geom_vline(xintercept = c(0, 1, 10, 20, 50))
    return(list(p1, p2))
}

## DLPFC lrt
gene_explore(ef_obj_dlpfc, de_lrt_dlpfc, get_mods(colData(rse[[1]]), only_pd = TRUE))

## ENSG00000227232.5 
##                 0

## [[1]]

## 
## [[2]]

## DLPFC odp
gene_explore(ef_obj_dlpfc, de_odp_dlpfc, get_mods(colData(rse[[1]]), only_pd = TRUE))

## [1] 7.822172e-07

## [[1]]

## 
## [[2]]

## HIPPO lrt
gene_explore(ef_obj_hippo, de_lrt_hippo, get_mods(colData(rse[[2]]), only_pd = TRUE))

## ENSG00000227232.5 
##                 0

## [[1]]

## 
## [[2]]

## HIPPO odp
gene_explore(ef_obj_hippo, de_odp_hippo, get_mods(colData(rse[[2]]), only_pd = TRUE))

## [1] 2.56492e-07

## [[1]]

## 
## [[2]]

## both
gene_explore(ef_obj, de_lrt, pd_all)

## ENSG00000217801.9 
##                 0

## [[1]]

## 
## [[2]]

#gene_explore(ef_obj, de_odp, pd_all)

Limma

This portion uses limma to fit the region by age interaction model and calculate a global F statistic.

library('limma')

## Warning: package 'limma' was built under R version 3.4.2

## 
## Attaching package: 'limma'

## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA

library('edgeR')

## Warning: package 'edgeR' was built under R version 3.4.2

library('qvalue')

ind <- data.frame(
    brnum = c(colData(rse[[1]])$BrNum, colData(rse[[2]])$BrNum),
    rnum = c(colData(rse[[1]])$RNum, colData(rse[[2]])$RNum),
    stringsAsFactors = FALSE
)

brnum <- ind$brnum[match(rownames(pd_all), ind$rnum)]

design <- mods$mod


counts <- do.call(cbind, lapply(rse, assay))[keepIdx, ]
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
v <- voom(dge, design, plot = TRUE)

## Page 49 (50 in the pdf) of limmaUsersGuide()
## Didn't run in my laptop (likely needs more mem)
if(FALSE) {
    corfit <- duplicateCorrelation(v$E, design, block=brnum)
    corfit$consensus
    fit <- lmFit(v, design, correlation = corfit$consensus)
} else {
    fit <- lmFit(v, design)
}

fit <- eBayes(fit)

colnames(design)[grep(':', colnames(design))]

## [1] "Age:RegionHIPPO"    "RegionHIPPO:fetal"  "RegionHIPPO:birth" 
## [4] "RegionHIPPO:infant" "RegionHIPPO:child"  "RegionHIPPO:teen"  
## [7] "RegionHIPPO:adult"

top <- topTable(fit, coef = grep(':', colnames(design)), n = nrow(counts))


qval <- qvalue(top$P.Value)
summary(qval)

## 
## Call:
## qvalue(p = top$P.Value)
## 
## pi0: 0.06969427  
## 
## Cumulative number of significant calls:
## 
##           <1e-04 <0.001 <0.01 <0.025 <0.05  <0.1    <1
## p-value    15061  15874 16968  17512 17984 18510 21015
## q-value    15913  17064 18727  19613 20538 21015 21015
## local FDR  14580  15312 16411  16925 17408 17961 21015

hist(qval)

summary(qval$qvalues)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.000e+00 0.000e+00 0.000e+00 4.029e-03 6.745e-05 6.965e-02

sapply(c(1e-04, 0.001, 0.01, 0.025, 0.05, 0.1, 1), function(x) { sum(top$adj.P.Val < x)})

## [1] 14958 15771 16854 17404 17855 18399 21015

# log2FC <- fit$coefficients[, 2]
# pvalue <- fit$p.value[, 2]
# qvalue <- qvalue(pvalue)$qvalues

Reproducibility

## Reproducibility information
print('Reproducibility information:')

## [1] "Reproducibility information:"

Sys.time()

## [1] "2017-10-23 14:30:02 EDT"

proc.time()

##     user   system  elapsed 
## 1117.413  169.320 1335.220

options(width = 120)
session_info()

## Session info ----------------------------------------------------------------------------------------------------------

##  setting  value                       
##  version  R version 3.4.1 (2017-06-30)
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/New_York            
##  date     2017-10-23

## Packages --------------------------------------------------------------------------------------------------------------

##  package              * version  date       source         
##  acepack                1.4.1    2016-10-29 CRAN (R 3.4.0) 
##  annotate               1.55.0   2017-05-04 Bioconductor   
##  AnnotationDbi          1.39.4   2017-10-12 Bioconductor   
##  assertthat             0.2.0    2017-04-11 cran (@0.2.0)  
##  backports              1.1.1    2017-09-25 CRAN (R 3.4.2) 
##  base                 * 3.4.1    2017-07-17 local          
##  base64enc              0.1-3    2015-07-28 cran (@0.1-3)  
##  bindr                  0.1      2016-11-13 CRAN (R 3.4.0) 
##  bindrcpp               0.2      2017-06-17 CRAN (R 3.4.0) 
##  Biobase              * 2.37.2   2017-05-05 Bioconductor   
##  BiocGenerics         * 0.23.3   2017-10-07 Bioconductor   
##  BiocParallel           1.11.11  2017-10-12 Bioconductor   
##  biomaRt                2.33.4   2017-08-02 Bioconductor   
##  Biostrings             2.45.4   2017-08-28 Bioconductor   
##  bit                    1.1-12   2014-04-09 CRAN (R 3.4.0) 
##  bit64                  0.9-7    2017-05-08 CRAN (R 3.4.0) 
##  bitops                 1.0-6    2013-08-17 cran (@1.0-6)  
##  blob                   1.1.0    2017-06-17 CRAN (R 3.4.0) 
##  BSgenome               1.45.3   2017-09-15 Bioconductor   
##  bumphunter             1.18.0   2017-10-05 Bioconductor   
##  checkmate              1.8.4    2017-09-25 CRAN (R 3.4.2) 
##  cluster                2.0.6    2017-03-10 CRAN (R 3.4.1) 
##  codetools              0.2-15   2016-10-05 CRAN (R 3.4.1) 
##  colorspace             1.3-2    2016-12-14 cran (@1.3-2)  
##  compiler               3.4.1    2017-07-17 local          
##  corpcor                1.6.9    2017-04-01 CRAN (R 3.4.0) 
##  data.table             1.10.4-2 2017-10-12 CRAN (R 3.4.2) 
##  datasets             * 3.4.1    2017-07-17 local          
##  DBI                    0.7      2017-06-18 CRAN (R 3.4.0) 
##  DelayedArray         * 0.3.21   2017-09-27 Bioconductor   
##  derfinder              1.11.10  2017-10-12 Bioconductor   
##  derfinderHelper        1.11.4   2017-10-12 Bioconductor   
##  devtools             * 1.13.3   2017-08-02 CRAN (R 3.4.1) 
##  digest                 0.6.12   2017-01-27 CRAN (R 3.4.0) 
##  doRNG                  1.6.6    2017-04-10 CRAN (R 3.4.0) 
##  downloader             0.4      2015-07-09 CRAN (R 3.4.0) 
##  dplyr                  0.7.4    2017-09-28 CRAN (R 3.4.2) 
##  edge                 * 2.9.0    2017-05-10 Bioconductor   
##  edgeR                * 3.19.8   2017-10-13 Bioconductor   
##  evaluate               0.10.1   2017-06-24 CRAN (R 3.4.1) 
##  foreach                1.4.3    2015-10-13 CRAN (R 3.4.0) 
##  foreign                0.8-69   2017-06-22 CRAN (R 3.4.1) 
##  Formula                1.2-2    2017-07-10 CRAN (R 3.4.1) 
##  genefilter             1.59.0   2017-05-04 Bioconductor   
##  GenomeInfoDb         * 1.13.5   2017-10-05 Bioconductor   
##  GenomeInfoDbData       0.99.1   2017-07-17 Bioconductor   
##  GenomicAlignments      1.13.6   2017-09-15 Bioconductor   
##  GenomicFeatures        1.29.13  2017-10-13 Bioconductor   
##  GenomicFiles           1.13.10  2017-07-17 Bioconductor   
##  GenomicRanges        * 1.29.15  2017-10-05 Bioconductor   
##  GEOquery               2.45.2   2017-09-26 Bioconductor   
##  ggplot2              * 2.2.1    2016-12-30 CRAN (R 3.4.0) 
##  glue                   1.1.1    2017-06-21 CRAN (R 3.4.1) 
##  graphics             * 3.4.1    2017-07-17 local          
##  grDevices            * 3.4.1    2017-07-17 local          
##  grid                   3.4.1    2017-07-17 local          
##  gridExtra              2.3      2017-09-09 CRAN (R 3.4.1) 
##  gtable                 0.2.0    2016-02-26 CRAN (R 3.4.0) 
##  Hmisc                  4.0-3    2017-05-02 CRAN (R 3.4.0) 
##  hms                    0.3      2016-11-22 cran (@0.3)    
##  htmlTable              1.9      2017-01-26 CRAN (R 3.4.0) 
##  htmltools              0.3.6    2017-04-28 CRAN (R 3.4.0) 
##  htmlwidgets            0.9      2017-07-10 CRAN (R 3.4.1) 
##  httr                   1.3.1    2017-08-20 CRAN (R 3.4.1) 
##  IRanges              * 2.11.19  2017-10-08 Bioconductor   
##  iterators              1.0.8    2015-10-13 CRAN (R 3.4.0) 
##  jackstraw              1.1.1    2017-09-04 CRAN (R 3.4.1) 
##  jsonlite               1.5      2017-06-01 CRAN (R 3.4.0) 
##  knitr                  1.17     2017-08-10 CRAN (R 3.4.1) 
##  labeling               0.3      2014-08-23 cran (@0.3)    
##  lattice                0.20-35  2017-03-25 CRAN (R 3.4.1) 
##  latticeExtra           0.6-28   2016-02-09 CRAN (R 3.4.0) 
##  lazyeval               0.2.0    2016-06-12 cran (@0.2.0)  
##  lfa                    1.7.0    2017-05-10 Bioconductor   
##  limma                * 3.33.14  2017-10-12 Bioconductor   
##  lme4                   1.1-14   2017-09-27 CRAN (R 3.4.2) 
##  locfit                 1.5-9.1  2013-04-20 CRAN (R 3.4.0) 
##  magrittr               1.5      2014-11-22 cran (@1.5)    
##  MASS                   7.3-47   2017-02-26 CRAN (R 3.4.1) 
##  Matrix                 1.2-11   2017-08-16 CRAN (R 3.4.1) 
##  matrixStats          * 0.52.2   2017-04-14 CRAN (R 3.4.0) 
##  memoise                1.1.0    2017-04-21 CRAN (R 3.4.0) 
##  methods              * 3.4.1    2017-07-17 local          
##  mgcv                   1.8-22   2017-09-19 CRAN (R 3.4.2) 
##  minqa                  1.2.4    2014-10-09 CRAN (R 3.4.0) 
##  munsell                0.4.3    2016-02-13 cran (@0.4.3)  
##  nlme                   3.1-131  2017-02-06 CRAN (R 3.4.1) 
##  nloptr                 1.0.4    2014-08-04 CRAN (R 3.4.0) 
##  nnet                   7.3-12   2016-02-02 CRAN (R 3.4.1) 
##  parallel             * 3.4.1    2017-07-17 local          
##  pkgconfig              2.0.1    2017-03-21 CRAN (R 3.4.0) 
##  pkgmaker               0.22     2014-05-14 CRAN (R 3.4.0) 
##  plyr                   1.8.4    2016-06-08 cran (@1.8.4)  
##  prettyunits            1.0.2    2015-07-13 CRAN (R 3.4.0) 
##  progress               1.1.2    2016-12-14 CRAN (R 3.4.0) 
##  purrr                  0.2.3    2017-08-02 CRAN (R 3.4.1) 
##  qvalue               * 2.9.0    2017-05-04 Bioconductor   
##  R6                     2.2.2    2017-06-17 CRAN (R 3.4.0) 
##  RColorBrewer           1.1-2    2014-12-07 cran (@1.1-2)  
##  Rcpp                   0.12.13  2017-09-28 CRAN (R 3.4.1) 
##  RCurl                  1.95-4.8 2016-03-01 cran (@1.95-4.)
##  readr                  1.1.1    2017-05-16 CRAN (R 3.4.0) 
##  recount              * 1.3.12   2017-10-13 Bioconductor   
##  registry               0.3      2015-07-08 CRAN (R 3.4.0) 
##  rentrez                1.1.0    2017-06-01 CRAN (R 3.4.0) 
##  reshape2             * 1.4.2    2016-10-22 cran (@1.4.2)  
##  rlang                  0.1.2    2017-08-09 CRAN (R 3.4.1) 
##  rmarkdown            * 1.6      2017-06-15 CRAN (R 3.4.0) 
##  rngtools               1.2.4    2014-03-06 CRAN (R 3.4.0) 
##  rpart                  4.1-11   2017-03-13 CRAN (R 3.4.1) 
##  rprojroot              1.2      2017-01-16 cran (@1.2)    
##  Rsamtools              1.29.1   2017-08-19 Bioconductor   
##  RSQLite                2.0      2017-06-19 CRAN (R 3.4.1) 
##  rtracklayer            1.37.3   2017-07-22 Bioconductor   
##  S4Vectors            * 0.15.14  2017-10-16 Bioconductor   
##  scales                 0.5.0    2017-08-24 CRAN (R 3.4.1) 
##  snm                    1.25.0   2017-05-10 Bioconductor   
##  splines                3.4.1    2017-07-17 local          
##  stats                * 3.4.1    2017-07-17 local          
##  stats4               * 3.4.1    2017-07-17 local          
##  stringi                1.1.5    2017-04-07 cran (@1.1.5)  
##  stringr                1.2.0    2017-02-18 cran (@1.2.0)  
##  SummarizedExperiment * 1.7.10   2017-09-29 Bioconductor   
##  survival               2.41-3   2017-04-04 CRAN (R 3.4.1) 
##  sva                    3.25.4   2017-06-20 Bioconductor   
##  tibble                 1.3.4    2017-08-22 CRAN (R 3.4.1) 
##  tidyr                  0.7.2    2017-10-16 CRAN (R 3.4.2) 
##  tools                  3.4.1    2017-07-17 local          
##  utils                * 3.4.1    2017-07-17 local          
##  VariantAnnotation      1.23.9   2017-10-12 Bioconductor   
##  withr                  2.0.0    2017-07-28 CRAN (R 3.4.1) 
##  XML                    3.98-1.9 2017-06-19 CRAN (R 3.4.1) 
##  xml2                   1.1.1    2017-01-24 cran (@1.1.1)  
##  xtable                 1.8-2    2016-02-05 cran (@1.8-2)  
##  XVector                0.17.1   2017-08-19 Bioconductor   
##  yaml                   2.1.14   2016-11-12 cran (@2.1.14) 
##  zlibbioc               1.23.0   2017-05-04 Bioconductor